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ON THE EFFICIENT USE OF PREDICTOR-CORRECTOR 
METHODS IN THE NUMERICAL SOLUTION OF 
DIFFERENTIAL EQUATIONS 


SUMMARY 

Basic interpolation formulas by Hermite can be used to generate large classes of 
correctors to be used in predictor-corrector processes for the numerical solution of 
differential equations. Conditions for the numerical stability of predictor-corrector pro- 
cesses are known [ref. 2, p. 20; 4, p. 197; 5, p. 218] . Only stable processes should be 
used in any solution requiring more than a few steps. 

Predictor-corrector methods are not self-starting; therefore, they are harder to 
program than starting processes. If programmed so that there are few wasted applica- 
tions of the corrector, however, then the predictor-corrector processes will have a 
significant computing time advantage over starting methods. 

The authors of this study are faculty members of the Mathematics Department of 
Vanderbilt University, Nashville, Tennessee, serving NASA under contract NAS8-2559. 


SECTION I. INTRODUCTION 

The primary purpose of this paper is to show how a predictor-corrector process 
can compete with a starting (self-starting) process. 

Section II gives the derivation of Hermite interpolation formulas and their error 
terms. Also, the predictor-corrector process for solving differential equations is out- 
lined. Numerical stability and propagated error are discussed. Some specific correctors 
are derived, and comments on the choice of a method are made. 

In section HI, the programming of predictor-corrector methods is discussed, 
analyzed, and illustrated. Comparisons with starting methods are discussed and charted. 
The usual programming of predictor-corrector methods is modified so as to save com- 
puting time, with no loss in numerical stability. 

Roundoff errors are not discussed, but at least four M guarding M digits are main- 
tained in all computations. That is, four digits beyond the accuracy of the method are 
carried along. 



SECTION II. PREDICTOR-CORRECTOR METHODS 


A. SOME BASIC FORMULAS 

We present formulas for approximating a function f (x) when certain values 
of the function and its derivative are known. The function is assumed to be analytic in 
the interval under consideration, and the points at which the function is evaluated are 
regularly spaced along the x-axis. A general method of deriving such formulas is ex- 
plained, and several specific formulas (with error terms included) are given [ref. 4, 
p. 128; 5, p. 209] . 

Let y(x) be the function to be approximated, and let xj = x + ih, i = 0, 1 

Denote the actual values of y(xj) and y'(xi) by zi and z J j. Let yi, y'j be approximations. 
We shall consider formulas of the type 


y n + l = 2 <Vj + “VJ). .<‘> 

1 = 0 


In the derivation of specific formulas of type (i) , the following is important: 

Theorem. If (1) is exact (that is, y . , = z , ) for y (x) = x m 
9 J n + 1 n + 1 J 7 

(m a positive integer) , then (1) is exact for any polynomial y(x) of degree not exceeding 

m. 


exact for y(x) = x 


Proof : First it is proved that if (1) is exact for y(x) = x 
m-1 

. Suppose, 


m 


then it is 


jx + (n + 1) b] m 
is an identity in x and h. 



Equating coefficients of x m a h a gives m + 


1 


equations 


2 


for a = 0, 


n 


i = Z A i 

i = 0 


n + i = Yj ( iA. + B. ) for a = i, and 
i = 0 ' 1 


(2) 


(n+ i) 


a v /m\ .a . , ( m — 1 A . a - 1 _ 

= >. i A. + , mi B. 

fJ \aj i \a> - i) i 


for a = 2, 3, . . . , m. 


Here | m ] is the usual binominal coefficient. Applying the rule 


a 


m 


three times to the last equation in (2) gives 


m fm-11 


a m-a \ a 


M (n+u “ -L [M‘“ a > + t-t) 


(m-l) i a_1 B 


for a = 2, 3, . . . , m-l. This implies (compare with (2)1 that ( 1 ) is exact for 

. . m-i ' J 

y(x) = x 


Next, observe that if formula (i) is exact for two functions f t and f 2 , then it is 
exact for ^ fj + C 2 f 2 , where Cj, C 2 are arbitrary constants. Hence the exactness of 

formula (1) for x m implies its exactness for any linear combination of x m , x m . . . 
x°, and the theorem is proved. 


A formula (1) is said to be of order m if it is exact for a polynomial of degree 
m but not exact for a polynomial of degree m + 1. 

In formula (1) the number of "back points" is n + i. The requirement that the 
2n “h 1 

formula be exact for z = x is met by choosing A^, Bj (i = 0, . . . , n) to satisfy a 

system of 2n + 2 linear equations. 
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In order to interpret formula (1) geometrically, suppose that the values of a 
function z and its derivative z' are known for x 0 , x 0 + h, . . . , x 0 + nh. If formula (1) is 
of order 2n + 1 , then the application has the effect of fitting a polynomial curve of degree 
2n + i through the n + 1 given points so that at each point the polynomial and the function 
z have the same derivative. Then y n + i is an approximation to z n + This is the well 
known "Hermite" interpolation [ref. 4, p. 96]. 

Next are listed some Hermite interpolation formulas. They have the form of 
formula (i) , with order 2n + i. In each case the number of back points is n + 1. The 
error terms, given in brackets, are discussed later in this section. 


yi = z 0 + h z 0 . 



(3) 


y 2 = 5z 0 - 4z t + h ( 2 zq + 4zi) . 



(4) 


y 3 = 10 z 0 + 9 zj - 18 z 2 + h (3 z{ + 18 z[ + 9z 2 ). 



(5) 


3y 4 = 47 z 0 + 192 z x - 108 z 2 - 128 z 3 + h (12 z 3 + 144 z{ + 216 z \ + 48 zj) . 


£ Z ( 8> 

70 



( 6 ) 


6y s = 131 z 0 + 1150 Zj + 600 z 2 - 1400 z 3 - 475 z 4 + h (30 zj+ 600 z{ 


+ 1800 z 2 + 1200 z 3 + 150 z\) ■ 


' h 10 
252 


z (lo) (s) 


(7) 


The derivations of the preceding formulas are straightforward. Formula (3) is 
recognized as a Taylor expansion. For the derivations of the other four, solve equations 
(2) for Ai, Bi. For example, take n = 1 and m = 3. Then equation (2) becomes 
Ag + A-l = i, Aj + Bg + B^ = 2, 3 + 6 B^ = 12, Aj + 3 B^ — 8. Hence Ag = 5, A^ = -4, 

B 0 = 2, Bi = 4, and formula (4) is obtained. Formulas (5) , (6) , (7) are derived 
similarly. 

The justification of the error terms is outlined in Hamming [ref. 4, p. 103] . The 
general error term is: 


4 



f (n + 1) ! 1 2 h 
( 2n + 2) ! 


2n + 2 


(s). 


( 8 ) 


(2n + 2) 
z 


The error terms in formulas (3) , . . . , (7) are found by taking n = 0, . . . , 4. These 
error terms represent the so-called truncation error and should not be confused with 
roundoff error. 

fe. THE NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS BY PREDIC- 
TOR-CORRECTOR METHODS 

The literature on the numerical solution of differential equations is vast 
( ref. 1; 4, p. 403; 5, p. 389] . Consider the differential equation y' = f(x, y) with initial 
condition y(x 0 ) = y 0 . There are various methods for producing points which lie on the 
solution curve [ref. i, p. 79; 4, p. 212; 9] . The Runge-Kutta method is widely used. 
Recently E. B. Shanks has developed similar processes which are markedly more 
efficient from the standpoint of accuracy and computing time. Shanks' formulas are 
self-starting, and investigations of the present writers and others show that these 
formulas are ideally suited for problems in which computing time, accuracy, and error 
propagation are important. 

Predictor-corrector methods have the disadvantage of not being self-starting. 
Before a predictor-corrector method requiring n back points can be used, a starting 
process must be applied n - 1 times. (All processes described in this paper use one of 
Shanks' formulas as a starting procedure. ) It is obvious that predictor-corrector pro- 
cesses are harder to program than starting processes. The only motivation in investiga- 
ting predictor-corrector processes is the possibility of gaining a time advantage in a 
particular problem or class of problems. The degree of this time advantage is discussed 
in sections which follow. 

Let (x 0 , y 0 ) , , (x m , y m ) be m + i regularly spaced starting points for the 

equation y' = f (x, y) . A formula 


Y n + i 


m 


Z 

i = 0 


(A, 


y ■ + 

•’n-i 


h B. 


y ' .) 

J n-r 


(9) 


can be used repeatedly to produce new points. Of course, at each point the value of the 
derivative must be computed by using y '. = f (x^, y^) . Usually it is wiser to use 

formula (9) each step to give a tentative or predicted value of y n + y' + ^ and then 

to correct the value of y^ + by repeated applications of a formula 
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+ 


(10) 


m m 

y . = Y. C. y . + h Y D. y' 

+ 1 H 1 J n-i f-* . l J n - l 

1 = 0 l = -1 


In this process, formula ( 9) is called the. "predictor" and formula (10) the "corrector. " 
In one step, formula (9) is used only once, but formula (10) may be used more than 
once. Note that for each application of the corrector, the function f(x,y) must be 
evaluated. The number of corrector applications per step is a significant factor in any 
time study. 

At first glance, one may be tempted to use only the predictor. This would cer- 
tainly give an advantage in computing time, but the disadvantage is well-known. Errors 
in yi are propagated unfavorably. This is extremely serious and usually prohibits the 
use of the predictor alone if the solution is to be extended for more than a few steps. 

For example, the size of the coefficients in the predictor (see formula (5)). 


y . = lOy + 9y , - 18y + h [ 3y ' n + I8y' . + 9y' 1 

•’n + 1 ,y n-2 J n-i J n ^ 3 n-2 J n-1 J nj 

indicates how unfavorably errors in y^ are propagated. 

So we use a predictor-corrector combination. Furthermore, the corrector must 
be selected with care, or else the unfavorable propagation of error will persist. 


C. NUMERICAL STABILITY IN PREDICTOR-CORRECTOR METHODS 

Let (x 0 , y 0 ) , . . . , (x m , y m ) be starting points for y' = f(x, y) , and apply 
the predictor-corrector pair 


m . . 

= / , I A. y . + h B. y ' .1 

n + 1 p l i J n-i l n - l I , 


( 11 ) 


m 


m 


Y . = Y C. y . + h y D. y' 

^n + 1 l J n-i H . l J n-i 

i = 0 l = -1 


(12) 


Suppose that in each step the corrector is applied repeatedly until two successive values 

of v . are equal. (Later this will be weakened, and it will be required only that two 
n + 1 

successive values of y R + ^ be "close" to each other.. A great deal of computing time 

will thus be saved. ) 

6 



Take z = z(x) as the exact solution of z T = f(x, z) , and let e ^ = z^ - be the 
difference between the true solution and the solution given by equations (ii) , (12). Sub- 
tracting equation (12) from 


m 

Z n + 1 = 2 


i = 0 


z . + 
n-i 


m 

z 


D. z 


-1 


» + 
n-i 


R, 


gives 


€ 


n + 


m 



m 

C. € . + h y D. e 4 . + R . 
i n-i H . l n-i 

l = -1 


(13) 


By the mean value theorem, there is a number t between y^ and z^ such that 


= <x r t)e i • 


(14) 


since e . = f(x^, z^) - f(x^, y^) . We follow the usual procedure [ref. 4, p. 197] in 

0f 

assuming that h— = K and R are constants. Then h e T . =Ke., and this with equation 

8 y . 11 

(13) gives the nonhomogeneous difference equation 


m 


X (C. + KD.) e . + R = 0 
H . i i n-i 

i = -1 


(15) 


in e with C . = -1. 
n-i -1 

If the roots \ 0 , . . . , A m of the ’'characteristic’ 1 equation 


m 


y (C. + K D.) 1 = 0 

H . l l 

i= -1 


( 16 ) 
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are distinct, then the general solution (e 


n +1’ ' ' ' ’ € n-m* ° f oc l uation < 15 ) is given by 


e 


k 




+ r. 


(17) 


where 


(r. 


n + r 


r n-m^ * s a P ar i-i cu ^ ar solution of equation (15) . 


From equation (17) it is seen that it is desirable to have the characteristic roots 
A 0 less than 1 in absolute value so that the propagated error will be favorable. 

(This last statement is not changed even if the characteristic roots are not all distinct. 
For a discussion see [ref. 5, p. 213] .) 

If h = 0, then K = 0 and one characteristic root is 1. For small h + 0, one 
characteristic root will be close to 1, but we require that the other characteristic roots 
be within the complex unit circle. It is impossible to predict in all situations what the 
range of values of K will be. Usually j k| is kept smaller than 0. 4 [ ref. 3, p. 7; 4, 
p. 198] . 

We will use only correctors such that, for h = 0, the characteristic roots (except 
1) are within the unit circle. There are many such formulas, but unless an investigation 
concerning the characteristic roots for h * 0 is made, Adams' formulas appear safer 
from the stability standpoint. In some instances, there are more efficient formulas. 

For h » 0 (h approximately equal to 0) , Fehlberg [ ref. 3, p. 5] has given a 
method for comparing the error propagation of two correctors. Since the authors have 
found this result of great use, it is reproduced here. Set h = 0 in equation (13) to get 


“n + 1 


m 

2 


C. e 
l n - 


+ C h 


m + 


3 (m 

y 


+ 3) 


(s) 


(18) 


and suppose there is a constant E defined by 


€ . = E k 
k 


. m + 3 
h 


(m + 3) 


(s) . 


(19) 


This is the assumption that the errors are almost in arithmetic progression. Now 
equations (18), (19) give 


8 



+ C, 


m 

E (n + 1) = Yj C. E (n - i) 
i = 0 1 


if y 


(m + 3) 


is assumed constant. 


m 

Since ^ 

i = 0 


C. 

l 


1 , we have 


E 



C. 

l 


( 20 ) 


If two correctors are used on the same differential equation, the long-run pro- 
pagated truncation errors will be proportional to the corresponding values of E, pro- 
vided the correctors are stable. The result is of great use in choosing a corrector. 


D. SOME SPECIFIC CORRECTORS 

Correctors of orders 3, 5, 7, 9 are presented in this section. Each 
corrector is a linear combination of formulas of types (4) through (7) . 

First we obtain correctors of order 3. Formula (4) gives: 

h 4 (4) 

Y2 = 5y 0 - 4yj + h (2y& + 4yJ) + y' ( s l) > ( 21 ) 


and reversing the order of the points (x 0 , y 0 ) , (x t , yj) , (x 2 , y 2 ) gives: 


Yo = 5y 2 - 4 yi - h ( 2y 2 + 4yJ) 



(s 2 ) 


( 22 ) 


A general formula for y 2 is obtained by multiplying equations (21) , (22) by 5b + 1, b, 
respectively, and then adding. The result is: 


9 



(23) 


y 2 = ( -24b -4) yi + (24b + 5) y 0 + h [ -2b yj + ( 16b + 4) y] + ( 10b + 2) yj] 

+ C h 4 y^ (s) , 

where C « b +4 , and s is between x 0 and x 2 . For h = 0, the characteristic roots of 
equation (23) are 1, -24b - 5. For numerical stability we require that < b < . 

Hence we have an infinite number of acceptable correctors. If we ignore the effect of 
roundoff error, the choice of a corrector (in this case, the choice of b) depends on two 
things: (1) the size of the step, or local, truncation error, Ch 4 y( 4 ) (s); (2) the pro- 
pagated truncation error. It is to be expected that there is no best formula for all 
problems. 

At one end, b = -— , of the stability interval, equation (23) becomes Simpson's 

rule. At b = - — , the result is 
4 

y 2 = 2y t - y 0 + h (1 y 2 y'o) . 


Simpson's rule has an attractive error term, - ^ h 5 y^ (s) , but is likely to be unstable 
[ref. 2, p. 37]. 

For b = , the midpoint of the stability interval, equation (23) becomes Adams f 

formula 

y 2 = yi + (5y' 2 + 8yi -yi) , (C = - £-) . (24) 

For b = - — , the harmonic mean of - -i and -_i , equation (23) becomes 
5 4. 6 

y 2 = f yi + |y 0 + | <2yb + 4y]) , (C = - ~) . (25) 


Computation of the ratio of E-values of equations (24) , (25) indicates that equation (25) 
will have the more favorable propagation of error. Numerical examples support this 
conclusion [ref. 3, p. 9] . 

We investigated one other formula from the infinite number of third-order 

formulas. Taking b = - — in equation (23) gives 

600 
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y2 = 25 yi + 25 y2 + 300 (l09y 2 + 328 y 'i + 55 y 0> 

(C » -0.015) . (26) 

The error propagated by formula (26) is about one-third of the error propagated by 
equation (25). This agrees with the ratio of the E-values. 

Similar formulas of orders 5, 7, and 9 are now given. For the derivations, see 
Appendix I. 


1 A x i 9 

Yi = 16 Y3 i Y2 4 Yl 16 y ° 

h (3703 yj + I5,518y' 3 + 6168 y^ + lO^Syi + 1873 y' 0 ) 

11,520 

+ C h 6 y (6) (s) , (27) 

where C » -0. 00725. For h = 0, the characteristic roots are approximately 1, -0. 800, 
-0. 069 ± 0. 836 i. The largest absolute value of a root is about 0. 84. A glance at the y- 
coefficients shows that they are nearly in geometric progression. 

A similar formula of order 7 is 


ye 


1 1 1 1 1 33 

M ys + 32 Y4 + 16 ys + 8 Y2 + 4 + 64 yo 


430, 080 


( 128, 627 y\ + 642, 168 yj + 130, 167 y| 


+ 693, 632 y 3 + 142, 137 y{ + 399, 240 y| + 61,469 yj) 


,8 ,.( 8 ) 


+ C h 8 y' 7 (S) , 


(28) 


where C » - 0. 00497. For h = 0, the characteristic roots have absolute values approxi- 
mately i, 0.85, 0.91, 0.91, 0.86, 0.86. 
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A formula of order 9 is: 


2, 560, 016 y 8 = 9784 y 7 + 20, 133 y 6 + 4i,040 y 5 

+ 79, 775 y 4 + 159, 816 y 3 + 319,691 y 2 + 639, 792 yi 

+ 1,289,985 y 0 + h (725, 340 y^ + 4, 150, 740 y} 

- 280, 710 y 8 + 6, 541, 620 yj - 1,808, 250 y| 

+ 5, 630, 940 y^ + 244, 290 y 2 + 2, 458, 620 y] 

+ 345, 330 y o ) + C h 10 y (10) (s) , (29) 

where C » - 0. 00361. The absolute values of the characteristic roots (for h = 0) are 
approximately!, 0.88, 0.80,0.80, 0.90, 0.90, 0.95, 0.95. 

Formulas (26) , (27) , (28) , (29) are examples of formulas with low truncation 
error terms and a reasonably stable behavior. These formulas, as well as Adams' 
formulas, are used to test and illustrate programming techniques in section III. Unless 
special circumstances arise, the authors prefer Adams' formulas because of their ex- 
cellent stability behavior [ ref. 8, p. 29J . 

E. CHOICE OF A METHOD 

In case a great deal of computing time is at stake, a well-chosen predictor- 
corrector process is likely to be better than a starting process alone. Of course, the 
former process is harder to program, but it can save time. If computing time is not 
critical, then starting methods, especially those developed by Shanks, are recommended. 
Also, it should be remarked that Shanks' formulas seem to be applicable to a large class 
of equations at large step sizes, without loss of stability. This is an advantage not shared 
by all predictor-corrector processes. 

If a predictor-corrector is used, the appropriate characteristic roots should be 
within the unit circle. If an estimate of K (see section II - C) is available, the charac- 
teristic roots should be examined further. If no such estimate is available, the appro- 
priate Adams formula is recommended. However, it will be shown in sections which 
follow that in many cases other correctors are better. 

The corrector, rather than the predictor, controls the error and stability 
behavior of these processes. A fifth-order predictor ( formula (5), section II) is used 
in nearly all of the processes reported in this paper. A few other predictors were tried, 
but they proved inferior to the one just mentioned. However, our experiments with 
different predictors were few, and it may well be that a small gain in computing efficiency 
would follow a more careful matching of predictor with corrector. 
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This report does not include a comparison with methods designed specifically for 
second (or higher) order differential equations. The processes presented are used to 
solve second-order equations in the traditional way, that is, by solving two first-order 
equations. It is emphasized that the final program presented can be used to solve up to 
twenty first-order equations, or up to ten second-order equations, etc. 


SECTION III. PROGRAMMING PREDICTOR-CORRECTOR METHODS 


A. INTRODUCTORY REMARKS 

Predictor-corrector methods usually require more than one back point. 
Consequently, a starter (starting process) is used to produce, from an initial point, 
the additional values needed for the predictor-corrector. 

Special mention is now made concerning programming methods. All programs 
were done by an IBM 7072. Fortran, with certain modifications, was used. It has been 
observed that roundoff error can affect the last four digits on runs of fewer than 2, 000 
steps. Hence, single precision (8-digit accuracy) does not allow a genuine comparison 
of numerical methods, and double precision (16 digits) is used. Double precision 
statements at the Vanderbilt computer are in modified subroutine form. In a Fortran 
program, these subroutines are called by the statements 

R = FMDF (Ai, A2, Bi , B2, Cl, C2) 

R = FADF (Al, A2, Bi, B2, Cl, C2) 

R = FDDF (Al, A2, Bi, B2, Cl, C2) 

These in turn produce the following: 


Cl 

+ 

C2 = (Al 

+ 

A2) 

(BI + B2) 

Cl 

+ 

C2 = (Al 

+ 

A2) 

+ (BI + B2) 

Cl 

+ 

C2 = (Al 

+ 

A2) 

/ (BI + B2) 


In each case Ai , Bi, Cl are the first eight digits and A2, B2, C2 the second eight digits 
of Al + A2, Bi + B2, Ci + C2, respectively. To illustrate, if Ai + A2 = 0. 12234 
55678 88999 9 * 10 15 , then Al = 0.12234 556 x 10 15 andA2 = 0 . 7 8 8 8 9 9 99 x 10 7 . 
Double precision programming with these subroutines is quite cumbersome. Each 
multiplication, addition, or division must be programmed as a statement of the types 
represented above. This means the programming of a number of predictor-corrector 
methods would be greatly delayed if each predictor-corrector were matched with a 
different starter. However, in any single process, it is desirable to match the starter 
with the predictor-corrector so that they have the same order. 
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In the studies herein presented, the following starter by Shanks [ref. 9, p. 7] 
was used. 


k 0 = h f r , 

kjt = h f(x r + f h ’ y r + | k o) , 

k 2 = hf[x r + |h, y r + | (k 0 + 3 ki)] , 

k 3 = hf[x r + |h, y r + ^- (4 k 0 + 6 k! + 8 k 2 )] , 

k 4 = hf[x r + ^h, y r + (17 k 0 + 12 k 4 + 16 k 2 - 9 k 3 )] , 

k 5 = hf[x r + |h, y r + Y^(llk 0 + 12 k t - 32 k 2 + 9 k 3 + 36 k 4 )] , 

k 6 = h f[x + h, y r + (-5 k 0 - 12 kj - 128 k 2 + 81 k 3 - 108 k 4 + 216 k 5 )] , 

y ,=y +t^t (U k 0 - 64k 2 + 81 k 3 + 81 k 5 + 11 k 6 ) . (30) 

r + l r xZk) 

This starter is of order six (the error is of order h 7 y ; (s) , where s is between x r and 
x r + i) . Because many of the predictor-corrector methods were of orders higher than 
six, additional accuracy was obtained for the initial values by using a smaller step size 
for the starting process than for the continuing process. For example, if five initial 
values at h = 0. 01 are needed for a seventh order predictor-corrector, these can be pro- 
duced as follows. From the given first point, produce eight new points at h = 0. 005. 
Denote the even numbered new points as (x 4 , y 4 ) , (x 2 , y 2 ) , (x 3 , y 3 ) , (x 4 , y 4 ) , and denote 



2 2 

/ y \ 

we have C 2“ e h 7 y v ' (s) as the error from (x., y.) to (x ., y .). For h = 0. 01 

/ v 11 1111. 

this error is approximately 1. 5625 Cy ! (s) x 10 -16 . Thus we have points of sufficient 

accuracy for a seventh order continuing process. Unfortunately, the technique does take 
additional time. 
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Another part of the program that is peculiar to our needs is that all coefficients in 
both predictor and corrector are in the form of data. Special branching is needed to 
generalize the usual programming. Such branching takes time. 


In the remainder of this section, we discuss the remaining sections and what we 
hope to accomplish. There are two basic ways of programming predictor-corrector 
methods. In one, the corrector is applied repeatedly at each step until two consecutive 
computed values are identical. In the. second basic technique, the corrector is applied a 
specific number of times per step. In section B, the problem of these two techniques are 
presented and illustrated. An appropriate convergence technique and variations thereof 
are presented in section C. The technique in section C saves considerable time over the 
technique of "complete convergence” in section B. Considerably more time can be 
saved, however, under the assumption that if convergence is approximate at one step, 
then it is approximate at another step. A program based on this assumption gave our 
best results. This, our final program, is presented in section D. 

For easy reference, we now list all computing processes and differential equations 
to be discussed in the remaining sections. The processes are called PI through P6, and 
the equations are called Ei through E4. In each process the predictor (see formula (5) ) 

y k+ 3 = - 18y k+ 2 + 9y k+l + 10y k + h (9y k+ 2 + 18y k + l + 3y k> 
is used. 

PI is the process in which Shanks' starting formula (30) is used as a continuing 
process. 

The other process are distinguished by their correctors, listed below. P3 and P5 
are by Adams. The source of the others is given in Appendix I. (The symbols y^ , yl 

are omitted in the right members of the following. The indices are in descending order. ) 


P2. 


k + 4 16 


= -- (1 + 2 + 4+ 9) + 


11, 520 


(3703 + 15, 518 + 6168 + 10, 898 + 1873) 


~ 167 

23, 040 


h 6 


( 6 ) 


(s). 


P3. 


y, ^*=(i + 0+ 0+ 0+ 0+0) +~ (19,087 + 65,112 - 46,461 + 37,504 

K + o oO, 4o0 

- 20, 211 + 6312 - 863) + h 8 y (8) (s) . 
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P4 - y k +6 = M <1 + 2 + 4+8 + 1 S+ 33 )+ iii^ <128,627 + 642,168 

+ 130,167 + 693,632 + 142,137 + 399,240 + 61,469) + h 8 y (8) (s) . 

P5. y k + 8 = ( 1 + 0+0+0+0+0+0+ °) + 3 628 800 < 1 ’ 111 ’ 267 + 4,137,094 
- 3,449,594 + 3,285,358 - 2,145,620 + 836,338 - 136,214 - 17,126 

+ 7297) + (-. 00936)h 9 (s). ( This is a modified Adams corrector. ) 

P6 ' y k + s ~ 2 560 016 < 9784 + 20 > 133 + 41 > 040 + 79 > 775 + l59 > 816 + 319 > 691 

•u 

+ 639,792 + 1,289,985) + ~ - (725,340 + 4,150,740 - 280,710 

2 , 5b0, 01 b 

+ 6,541,620 - 1,808,250 + 5,630,940 + 244,290 + 2,458,620+ 345,330) 

+ (- 0 . 00361) h 10 y^°^ (s). 


The equations to be solved numerically are given next. Included are the initial 
values, the interval of the solution, and the correct solution. 

El. y’ = y. Initial value y(0) = 1. Interval of solution x = 0 to x = 18. Correct 
solution y = e x , y( 18) = 65659 969.13733 05111 38786. 


E2. y"- = - X ^ X y ^ f • Initial values y(i) = y T - (1) =1. Interval x = 1 to x = 19. 

Correct solution y = y\J i + 2 log x , y' = ~ , y(19) = 2.62466 72090 63443, 


y' (19) = 0. 02005 26675 40335 10. 

E3. y' = -y. Initial value y(0) = 1. Interval x = 0 to x = 18. Correct solution 
y = e" X , y( 18) = 0.15229 97974 47126 284 xlO -7 ! 
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E4. y' = -2 x y 2 . Initial value y(G) = i. Interval x = 0 to x = 18. Correct 
solution y = (1 + x 2 )" 1 , y(!8) = 0. 00307 69230 76923 07692 31. 


Future reference to the preceding is made according to the labels. For example, 
P2 - Ei means the process P2 used to solve equation El. 


The equation E2 is solved by solving the pair z T 


XZ + y 
(xy)2 


z. 


B. THE USUAL PREDICTOR-CORRECTOR TECHNIQUES 


Both of the two basic techniques for programming predictor-corrector 
methods have definite problems. In solving differential equations, our goal is always to 
obtain a given accuracy with the smallest computing time. Because elegant starting pro- 
cesses exist, it becomes essential to find an efficient and dependable predictor-corrector 
program to compete effectively with these starting methods. Neither of the common 
methods is both efficient and dependable. 



FIGURE i. USING A SPECIFIC NUMBER OF CORRECTIONS 


Figure i is an abbreviated flow chart for one step of a predictor-corrector 
method. (LL is the number of corrections desired. ) Such a technique can be quite 
efficient as far as time alone is concerned. However, the method lacks predictability. 

If relative error builds up arithmetically (normally it appears to do so) , it is not apparent 
what this error should be. At any given step, actual convergence may not occur until 
the tenth correction or later. Error approximations that arise from Taylor’s series are 
based on actual convergence. As Table I shows, errors may change a great deal as the 
number (LL) of applications of the corrector per step is varied. 
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TABLE i. ERRORS FOR DIFFERENT NUMBERS OF CORRECTIONS. P4 - El. 
STARTERS PRODUCED AT ONE-HALF h 


Number 

of 

Corrections 


h 



0. 08 

0. 20 

1 

0. 1498 x 10-3 

0. 2976 x 10 _1 

2 

0.4263 x 10" 6 

0. 2252 x 10~ 3 
0. 2592 x 10 -2 

3 

0.4462 x 10 -5 


The unpredictability of this technique makes it impractical in general. 



FIGURE 2. CONVERGENCE AT EACH STEP 


The abbreviated flow chart of Figure 2 illustrates another technique. As we have 
mentioned, the error terms which are calculated from Taylor's series are valid. Hence, 
this technique is predictable as to error. For this reason, this technique is used very 
often. The drawback is its lack of efficiency. As shown in Table 2, considerable time 
is used to "converge to zero" in many cases. Hence, this method is impractical. 
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Note: By "complete convergence” or "convergence to zero" we mean application 
of the corrector each step until two successive values of are identical. 


TABLE 2. P4-E1. STARTERS PRODUCED AT ONE -HALF h 



4 Corrections Per Step 

Complete Convergence 

h 

Error 

Time (min) 

Error 

Time (min) 

0. 15 

2. 015 

0. 2920 

2. 019 

0. 4072 

0. 18 

7. 196 

0. 2477 

7. 215 

0. 3583 

0. 20 

14. 99 

0. 2261 

15. 04 

0. 3345 

0. 24 

53. 11 

0. 1931 

53. 36 

0. 2964 

0. 30 

246. 9 

0. 1603 

248. 8 

0. 2561 


C. APPROXIMATE ERROR RATIO CONVERGENCE 

Consider the seventh order process P4. In this, the error term E ^ is 

/ o \ 

less than or equal to (285/57, 344) h 8 y ; (s) , where s lies between x. and x 
/ 8 \ ( 8 ) 1 1 
Normally, y w (s) is not known because y (x) is seldom known and s is not predictible. 

However, it may be possible to approximate + ^ with some value + ^ calculated 

from what is known. Let (y. . ) . and (y. , , ) . , , be two successive approximations 

^l+ij + i j + i 

of y. + ^ with a given corrector. We make the following definition of variations on con- 
vergence. 


error if 


Definition i. (y. , . ) converges to (y. , , ) . within a ratio r of the 
— J i+t n J i+i j 


(y i + 1> j - ( yi + i>j + i 


r E i + i 


Definition 2 . (y. + ) converges to (y. + . within a ratio r of the 

i in i I ^ J . 

approximate error if 


<y i+ i>j ■ <y i+ 1*) + 1 


rE Li 
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Let a given predictor-corrector technique require n back points. Define, for 


1 = n ’ y i+l 
E. 


= (y. , , ) ., where (y. , , ) 
J i + i j J i + I n 


(y. , ^ ) • within a ratio r of the error 
i + 1 J 


. . By using this approximation of y in place of y , we can calculate y. 

Continuing the process inductively, we can approximate the solution given by the technique 
requiring complete convergence at each step. This new process will be called "error 
ratio convergence. ” If, instead of Definition 1, we use Definition 2 in approximating 


y,- 


i + 1’ 


we will call the process of solving the equation "approximate error ratio conver- 


gence. " This process has the following flow chart for one step. 


j^ y i + l^j “ ^ y i + l^j + i 


rE Uf 



FIGURE 3. APPROXIMATE ERROR RATIO CONVERGENCE 


We now give a useful way of obtaining E'. For the present, assume that the 
y i’ y i are exact ‘ For y = y< x )> define backward differences Vy. = y^ - y. ^ , and 

V j * y^ = V^y^ - V^y^ ^ , in the usual manner. Use h n V n y^ as an approximation for 

y^(x^) [ref. 6, p. 128]. Normally, an n-th order process requires n-1 back points 

and has an error E. of the form 
i 


E. 

l 


„ ,n + 1 (n+1). . 

Ch y (sj) forsjinfx. _ n + 1 , x.). 
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We can then define E^' by: 

E! - C h V n y \ for i — n, 

1 i 

where the exact values yl^ are replaced by the approximate values y! , which are calculated 
by the process being used. This definition of E' is quite useful. 

For an n-th order process, the calculation of E' requires n back points, one more 

than would be needed without this calculation. In all of the following sections E! is cal- 
• 1 

culated by using y! for j < i, and the predicted value of y^ in place of y! . 

We will use E^ in three ways. First E| + ^ can be calculated at each step with 
y. + i defined as (y i + ^ , where (y. + 1 > n — (y. + ^ within a ratio r of E'. This 
calculation of E! at each step does consume needless time, as we will demonstrate. In 
Table 3 we see the result of using one process with E! calculated at each step, (r is 
taken to be 0. 04. ) 


TABLE 3. P4 - El. STARTERS PRODUCED AT ONE-HALF h 



Approximate Error Ratio 
Convergence r = 0. 04 

Complete Convergence 
Each Step 

h 

Error 

Time (min) 

Error 

Time (min) 

0. 12 
0. 15 

0.4232 

0. 4989 

0.4232 

0. 4770 

2. 018 

0. 4039 

2. 019 

0. 4072 

0. 18 

7. 214 

0. 3397 

7. 215 

0. 3583 

0. 20 

15. 03 

0. 3083 

15. 04 

0. 3345 

0. 24 

53. 34 

0. 2600 

53. 36 

0. 2964 

0. 30 

248. 6 

0. 2125 

248. 8 

0. 2561 


It is clear that the calculation of Ej does require some time. Also, an additional correc- 
tion is used needlessly to show that (y^ + ^) n — (y. + within a ratio r of E| + For 
example, if (y. + ^ — (y. + 1 > 5 within a ratio r of E| + ^ we must find (y. + 1 ) 6 to 
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show that Definition (2) is satisfied. Normally, if (y^ + (y^ + 5 within a ratio r of 

E { + 1 ’ then (y 4 + 1 ) n -' (y i+ 1 ) e within a ratio rofE| + ^ 

It is possible to reduce this needless waste of time by further changes in the pro- 
gram. One change is given in this section, another in the next section. 

Define F., F! by 
11 J 

F. = (y.E^/yi , 

F! = (y.Ei)/ yi . 

If we use F! instead of E! in the technique just described, additional time is saved. 

Essentially, the use of F^ means that we are taking the relative error at the first step as 

an approximation for the relative error at the i-th step. The use of F! is similarly inter- 
preted. 

The following Table 4 shows the time gained by using Fh It is emphasized that 
the values of y^ produced by the two programs are identical. 


TABLE 4. P4-E1. STARTERS PRODUCED AT ONE-HALF h. r = 0 . 04 



Computing Time in Min 

h 

Using E'j 

Using F! 

0 . 12 

0.4989 

0. 4058 

0. 15 

0.4039 

0. 3303 

0 . 18 

0. 3397 

0. 2794 

0 . 20 

0. 3083 

0. 2541 

0. 24 

0. 2600 

0 . 2161 

0. 30 

0. 2125 

0.1784 | 


As before, it is still necessary to correct each step an additional time in order 
to make the decision that enough corrections have been made. This will be remedied in 
the next section where the technique is modified again. In the new program, errors differ 
only slightly from those obtained with the preceding technique. 
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D. THE FINAL PROGRAM— AN ASSUMED APPROXIMATE ERROR RATIO 
CONVERGENCE 

The final program involves the additional, though trustworthy, assumption 
that if ( Yj ) n (Yj) m within a ratio r of the approximate error Ej, then, for all i, 

(y ± ) n (y.) m within a ratio r of the approximate error E! or else the deviation is so 

slight as to be insignificant. 

This means that (y^) m + ^ needs to be computed only for the first step in which 

the predictor-corrector is used. At the first step, the number of corrections is computed. 
From that point on, the number of corrections per step is constant. This gives our final 
program, an ''assumed approximate error ratio convergence" technique. Its flow chart 
is shown in figqre 4. 

In tables 5 and 6 we see not only the time saved, but also the change in error, 
when the final program is used in place of the one which requires complete convergence 
at each step. Comparison of table 6 with table 4 shows the time advantage of the final 
program over the two earlier methods. In tables 5 and 6, the "percent change in error" 
entries should not be misinterpreted. For example, the reading 2. 17 in table 5 comes 
from comparing two relative errors, each of which is approximately 1CT 6 . 

In one instance, the time for the complete convergence process exceeded fourteen 
minutes, while the time for the final program was 0. 5672 minutes. This however, was 
unusual. Tables 5 and 6 represent typical time advantages. 


23 




FIGURE 4. THE FINAL PROGRAM 
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TABLE 5. P2-E1. r = 0. 08. 



Complete Convergence 

Final Program 

Comparison 

h 

Error 

Time (min) 

Error 

Time (min) 

Percent 
Change in 
Error 

Percent 

Time 

Saved 

0. 12 

67. 65 

0.4331 

66. 18 

0. 2464 

2. 17 

43. 11 

0. 15 

208. 2 

0. 3703 

201. 3 

0. 2003 

3. 31 

45. 91 

0. 18 

521. 8 

0. 3222 

497. 6 

0. 1692 

4. 64 

47. 49 

0. 20 

887. 1 

0. 3033 

837. 3 

0. 1539 

5. 61 

49. 26 

0. 24 

2220 

, 0. 2678 

2048 

0. 1305 

7. 75 

51. 27 

0. 30 

■i 

6805 

0. 2314 

6726 

0. 1180 

1. 16 

49. 01 


TABLE 6. P4-E1. r = 0. 04. 



Complete Convergence 

Final Program 

Comparison 

h 

Error 

Time (min) 

Error 

Time (min) 

Percent 
Change in 
Error 

Percent 

Time 

Saved 

0. 12 

0.4232 

0.4770 

0. 4227 

0. 3570 

0. 11 

25. 16 

0. 15 

2. 019 

0.4072 

2. 015 

0. 2920 

0. 19 

28. 29 

0. 18 

7. 215 

0. 3583 

7. 196 

0. 2477 

0. 26 

30. 87 

0. 20 

15. 04 

0. 3345 

14. 99 

0. 2261 

0. 33 

32. 41 

0. 24 

53. 36 

0. 2964 

53. 11 

0. 1931 

0. 46 

34. 85 

0. 30 

248. 8 

0. 2561 

246. 9 

0. 1603 

0. 76 

37.41 


25 



E. A NOTE ON MATCHING STARTER WITH PREDICTOR-CORRECTOR 


One of Shanks' sixth order starters, equation (30), was used with every process 
reported in this paper. Whenever the corrector was of higher order than the starter, the 
step size for the starter was taken as a fraction of that of the continuing process. For 
correctors of order seven, the starters were produced at half the regular step size ( see 
section III - A) . Similarly, starters for a ninth-order corrector were produced at a fifth 
of the regular step size. This was an obvious waste of time for large step sizes. 

/ 

The error of the starter is C* h 7 y ; (s*) , while the error of an n-th order 

corrector is C 2 h 11 + ^ y + ^ (s 2 ) . For h = 0. 01 and n = 7 it is sufficient to compute 
the starting points at half the regular step size. The results for a ninth-order corrector 
at step sizes 0. 01 and 0. 02 are confused by roundoff error, even though double precision 

is used. For h = 0. 03 and n = 9, corrector error is 5.9049 (i0~ 16 ) C 2 y (s 2 ). For 

starters produced at h/5, the error is approximately 5(2. 79936) (10 -16 ) C t y ; (s^ or 

1. 39968 (10~ 15 ) Cj y^ 7 ( Sj) . The use of an eighth-order starter would reduce calculation 

time considerably. 

For larger step sizes the starters may be produced more quickly. If h = 0. 3 and 
n = 9, the corrector error is 5. 9049 (10 -6 ) C 2 y^ 10 ^ (s 2 ). If the starters are produced at 
h/2, then the starter error is approximately 2( 1. 7) (10 -6 ) Cj y 7 (sj). 

It may be that using the starting process for so many steps is unnecessary in 
some cases, as the following table indicates. This particular aspect deserves attention 
in future studies. The question, "How do inaccuracies in initial values affect the final 
results?" would be pertinent to such a study. 

The approximate values in Table 8 are useful in estimating the time wasted in 
producing initial values at a fraction of the regular step size. 
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TABLE 7. P5 - E3 



Starters Produced at 1/5 h 

Starters Produced at h 

h 

Error 

Time (min) 

Error 

Time (min) 

0. 03 

0. 1606 x i0' 19 

1. 2011 

0. 1325 x 10" 19 

1. 1164 

0. 06 

0. 5402 x 10 -18 

0. 6514 

0. 3014 x 10“ 18 

0. 5667 

0. 09 

0. 1478 x 10" 16 

0.4683 

0. 1058 x 10 -16 

0. 3834 

0. 12 

0. 1791 x 10" 15 

0. 4037 

0. 1467 x 10 -15 

0. 3189 

0. 18 

0. 5315 x 10 -14 

0. 3028 

0. 4729 x 10~ 14 

0. 2181 

0. 24 

0. 6150 x i0 -13 

0. 2520 

0. 5685 x 10 -13 

0. 1675 


TABLE 8. TIME USED-FOR PRODUCING STARTERS AT A FRACTION OF THE REGU- 
LAR STEP SIZE 


Equation 

Time (min) for One 
Starter Application 

Time Wasted in Program 

Seventh-order. 
Starters Pro- 
duced at 1/2 h 

Ninth-order. 
Starters Pro- 
duced at h/ 5 

El 

0. 002666$ 

0. 0160 

0. 0853 

E2 

0. 006342 

0. 0381 

0. 2029 

E3 

0. 0026743 

0.0160 

0. 0856 

E4 

0. 0035717 

0. 0214 

0. 1143 
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SECTION IV. CONCLUSIONS 


Our results indicate that in a large number of cases predictor-corrector processes 
may compete successfully with starting methods. However, unless care is taken to re- 
duce the number of wasted applications of the corrector, then the starting process alone 
is more efficient. 

The final program given in section HI - D demonstrates clearly that predictor- 
corrector processes can save significant time when compared with very efficient starting 
processes. 

There is no single corrector (of a given order) which is the best one for all 
differential equations. Two correctors may compete differently when applied to distinct 
equations or when used at different step sizes. The process with the smaller error term 
may be efficient and accurate in one situation, but unstable in another. 

If the number K of section II - C can be approximated in advance, it is possible to 
design a corrector specifically for a given problem. 
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APPENDIX I. DERIVATION OF FORMULAS 
The derivations of formulas (27) , (28) , and (29) are outlined: 

Equation (27) . From section I, equation (5) , 

^6 /g\ 

Y 3 = 10 y 0 + 9 yi - 18 y 2 + h (3 y{, + 18 y\ +9 yp + -j^ y (si) , 

h 6 ( 6) 

yo = io Y3 + 9 y 2 - 18 yi - h ( 3 yj + 18 yJ + 9 yl) + ^ y ( s 2 ) > 

^6 / g\ 

y 4 = 10 y! + 9 y 2 - 18 y 3 + h (3 yl + i 8 yj + 9 yj) + ^ y' ’ (s 3 ) > 

h 6 (6) 

Yi = 10 Y4 + 9 Y 3 - 18 Y 2 - h (3 yi + 18 y$ + 9 yj) + ^ y (s 4 ). 

Multiply the preceding equations in turn by 1873, -710, -2470, -3703, and then add. 
This leads to (27). 

Equa tion (28) . From section I, equation (6), 

3 y 4 = 47 y„ + 192 y 4 - 108 y 2 - 128 y 3 + h ( 12 yj + 144 y { + 216 y\ + 48 yj) , 


3 y 0 = 47 y 4 + 192 y 3 - . . . -h(12yi+... • ), 

3 ys = 47 y 4 + 192 y 2 - . . . +h(12y{+... ), 

3 yi = 47 y 5 + 192 y 4 - . . . -h(12y£+... ), 

3 y 6 = 47 y 2 + 192 y 3 - . . . +h(i2yj>+... ), 

3 y 2 = 47y 6 + 192 y 5 - ... -h(12yfc+... ). 
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Multiplying the preceding equations respectively by 184,407; 227,923; -103,472; -833,968; 
-884, 509; -385, 881 and then adding gives equation (28). 

Equation (29). From section I, equation (7), we have: 


6 y 5 = 131 y 0 + 1150 y t + 600 y 2 - 140.0 y 3 - 475 y 4 
+ h (30 y[> + 600 y| + 1800 y £ + 1200 y$ + 150 y4) , 


and similar equations for 6 yo, 6 y 6 , 6 y 4 , 6 y 7 , 6 y 2 , 6 y 3 , 6 y 3 . As in the derivation of 
(27), (28), use multipliers 11, 511; 36,326; 33,364; -47,718; -135,347; -160,883; 
-101,217; -24,178. 

It is noted that formulas (27) , (28) , (29) are linear combinations of sets of 
formulas obtained from section II - A. Of course, an infinite number of correctors can 
be so generated. The ones just presented have their y-coefficients in approximate geo- 
metric progression. This makes the error term attractive, and for h = 0 the character- 
istic roots are within the unit circle. At large step sizes, Adams' formulas are some- 
times preferred because they are less likely to become unstable. 

Adams' formulas may be derived also by taking linear combinations of formulas 
obtained from section II - A. Also, they are found in the literature [refs. 1, pp. 10, 536; 
5, pp. 194, 204]. They are given by 


'r + 


ry r + h Z B i vl y' 

1 1 ,• — n 


r + 1 


+ h 


p + 2 


B 


(P + 2) 


/ c \ 


The first few coefficients are B 0 =i, B 4 = - B 2 = - Bs = _ 24 ’ B4 = “ 7^0 ’ 

3 _ 863 _ 275 _ 33, 953 73, 647 

Bs " " 160 5 Be “ " 60,480 ’ B? “ " 24, 192' Bs " "3, 628, 800 ’ 9 9,331,200’ 


Equivalent formulas (without remainders listed) are: 


yi = yo + 2 (yt> + yi) > 

y2 = yi + ^ - yo + 8 yl + 5 y2) » 
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y3 = y2 + -^ (yo - 5 y{ + 19 y 2 + 9 y 3 ) , 



1 

l 

Y4 = y 3 + ( - 19 yo + i06 y{ - 264 y' 2 + 646 y 3 + 25i y^) , 

Y5 = y 4 + (27 yi, - 173 y[ + 482 y\ - 798 y 3 + 1427 yi - 1475 y' 5 ) , 

ye = y5+ 6Mt80 (_ 863 y 0 + 6312 Yi - 20 ’ 211 Y2 + 37, 504 y^ - 46,461 yi 
+ 65,112 y' 5 + 19,087 y ' e ) . 

A similar formula of order eight is listed as P5 in section III-A. 
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APPENDIX II. TABULATED RESULTS. 


In these tables we include the estimated loss of time mentioned in section III - E. 
Notations identifying processes and differential equations follow that of the text. For P2 
through P6, the final program was used (see section III) . The constant r is the one 
defined in section ILL 


TABLE A-l. RESULTS OF Ei. 


Process 

h 

Error 

Time (min) 

Starters 
Produced At 

Mins. 
Time Loss 


CO 

o 

0. 5636 x 10" 3 

1. 6075 




. 06 

0. 3443 x lO -1 

0. 8069 



Pi 

. 10 

0. 7100 x 10° 

0.4881 




. 12 

0. 2080 x 10 1 

0.4073 




. 18 

0. 2236 x 10 2 

0. 2741 




. 24 

0. 1186 x 10 s 

0. 2075 




. 03 

0. 5482 x 10" 1 

0. 8305 

h 

0 


. 06 

0. 2057 x 10 1 

0.4789 

If 

0 

P2 

. 10 

0. 2660 x 10 2 

0. 2936 

n 

0 

(r = . 08) 

. 12 

0. 6618 x 10 2 

0. 2466 

n 

0 


. 18 

0.4976 x 10 s 

0. 1692 

n 

0 


. 24 

0. 2048 x 10 4 

0. 1305 

ff 

0 


O 

CO 

0. 1948 x lO" 3 

1. 0475 

h/2 

0. 0160 

P3 

. 06 

0. 3231 x 10 _1 

0. 5953 

I? 

ff 

o 

11 

u 

. 09 

0. 5128 x 10° 

0.4070 

ff 

ff 
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TABLE Ai. RESULTS OF Ei. (Concluded) 


Error 


Starters Mins. 

Time (min) Produced At Time Loss 


. 12 

0. 3565 x 10 1 

00 

0. 5245 x 10 2 

. 24 

0. 3385 x 10 3 

00 

o 

0. 1219 x iO -4 

. 06 

0.3269 x 10“ 2 

. 10 

0. 1178 x 10° 

. 12 

0.4227 x 10° 

. 18 

0. 7196 x 10 1 

. 24 

0. 5311 x 10 2 

. 03 

0. 7016 x 10" 4 

. 06 

0. 1566 x 10 -2 

. 09 

0. 3710 x 10 -1 

. 12 

0. 3900 x 10° 

. 18 

0. 8225 x 10 1 

. 24 

0. 6763 x 10 2 

. 03 

0. 1540 x 10- 4 

. 06 

0. 3689 x 10- 4 

. 10 

0. 5471 x 10 -3 

. 12 

0. 2833 x 10 -2 

. 18 

0. 1059 x 10° 

. 24 

0. 1345 x 10 1 


0. 3128 
0. 2186 
0. 1720 
1. 2289 
0. 6850 
0.4234 
0. 3570 
0. 2477 
0. 1931 
1. 1992 
0. 6503 
0.4672 
0. 4028 
0. 3019 
0. 2513 
1. 5636 
0. 7744 
0. 5377 
0.4633 
0. 3414 
0. 2803 


0. 0160 


0. 0160 


0. 0853 


0. 0853 


TABLE A2. RESULTS OF E2 


Process 

h 

Error 

Time (min) 

Starters 
Produced At 

Time Loss 
(min) 


. 02 

0. 1845 x io- 9 

5. 7336 




. 04 

0. 1172 x 10 -7 

2. 8706 



Pi 

. 08 

0. 7258 x IO -6 

1. 4403 




. 12 

0. 7780 x 10" 5 

0. 9638 




. 18 

0. 7795 x IO -4 

0. 6467 




. 24 


0.4881 




. 03 

I 

2. 2150 

h 

0 


. 06 


1. 1203 

M 

0 

P2 

. 10 

0. 2685 x IO -4 

0. 6836 

rr 

0 

(r = . 08) 

. 12 

0. 2538 x IO -4 

0. 5733 

TT 

0 


. 18 

0. 9650 x 10 -4 

0. 3911 

" 

0 


. 24 

0. 5559 x IO -3 

0. 3000 

IT 

0 


. 03 

0. 6210 x 10" 6 

2. 2022 

h/2 

0.0381 


. 06 


1. 1361 

I! 

TT 

P3 

. 09 

0. 9929 x 10 _4 

0. 7808 

If 

TT 

(r = . 04) 

. 12 

0. 2725 x IO -3 

0. 6033 

ft 

TT 


. 18 

0. 8876 x IO -3 

0.4258 

TT 

TT 


. 24 

0. 1729 x 10“ 2 


IT 

TT 
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TABLE A2. RESULTS OF E 2 (Concluded) 


Process 

h 

Error 

Time (min) 

Starters 
Produced At 

Time Loss 
(min) 


. 03 

0. 3518 x 10“ 7 

2. 5664 

h/2 

0. 0381 


. 06 

0. 6126 x 10" 6 

1. 3164 

M 

ft 

P4 

1 o 

0.7187 x i0" 6 

0. 8175 

11 

IT 

o 

II 

. 12 

0. 1971 x 10 -5 

0. 6917 

! T 

M 


. 18 

0. 3483 x 10 -4 

0.4836 

1 f 

IT 


.24 

0. 1166 x 10 -3 

0. 3794 

M 

IT 


. 03 

0. 1472 x 10 -6 

2. 8772 

h/5 

0. 2029 


. 06 

0. 5140 x 10 -5 

1. 3867 

1! 

IT 

P5 

. 09 

0.2725 x 10 -4 

1. 0042 

IT 

II 

(r = . 04) 

. 12 

0. 7352 x 10" 4 

0. 8127 

If 

n 


. 18 

0. 2444 x 10“ 3 

0. 5683 

II 

IT 


.24 

0.4357 x 10~ 3 

0.4872 

IT 

IT 


. 04 

0. 9493 x 10 -8 

2. 5948 

h/5 

0. 2029 


. 06 

0. 3740 x 10~ 7 

1.6367 

IT 

IT 

P6 

. 10 

0. 6723 x 10~ 6 

1.0761 

II 

IT 

(r = . 04) 

.12 

0. 2375 x 10- 5 

0. 9344 

11 

IT 


. 18 

0. I899x 10~ 4 

0. 6472 

rr 

IT 


. 24 

0. 5845 x 10~ 4 

0. 5834 

it 

IT 
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TABLE A3. RESULTS OF E3 


Error 

0. 1382 x IQ" 18 


Starters Time Loss 
Time (min) Computed At (min) 

1. 6106 


0. 8964 x 10” 17 


0. 8083 


0. 1051 x 10 -15 
0. 6077 x 10“ 15 


0. 5414 


0.4078 


0. 7334 x 10' 14 


0. 2744 


0. 4366 x 10 -13 


0. 2075 


0. 1662 x 10" 16 


0. 8322 


0.4399 x 10 -15 


0.4800 


0. 3222 x 10" 14 
0. 1303 x 10 -13 


0. 3247 


0. 2472 


0. 2818 x 10“ 13 


0. 1698 


0. 2435 x 10 -9 
0. 5743 x 10 -19 
0. 9824 x 10 -17 



0. 1306 


1. 0488 


0. 5403 


0. 3709 


0 . 2861 


0. 2014 


0. 2350 x 10 









TABLE A3. RESULTS OF E3 (Concluded) 


Process 

h 

Error 

Time (min) 

Starters 
Computed At 

Time Loss 
(min) 


. 03 

0. 3890 x IO -20 

1. 2317 

h/2 

0. 0160 


. 06 

0.7173 x 10 -18 

0. 6309 

?t 

ft 

P4 

. 09 

0. 1194 x 10 -16 

0. 4303 

ft 

ft 

( r = . 04) 

. 12 

0. 1430 x 10 -15 

0. 3305 

n 

tf 


. 18 

0. 1296 x 10“ 10 

0. 2302 

I! 

ft 


. 24 

0.4631 x 10“ 8 

0. 1803 

ft 

ft 


. 03 

0. 1606 x 10 -19 

1. 2011 

h/5 

0. 0856 


. 06 

0. 5402 x 10 -18 

0. 6514 

It 

ft 

P5 

. 09 

0. 1478 x 10“ 16 

0.4683 

ft 

1 1 

(r = . 04) 

. 12 

0. 1791 x 10" i5 

0.4037 

ft 

ft 


. 18 

0. 5315 x 10 -14 

0. 3028 

If 

ft 


. 24 

0. 6150 x 10“ 13 

0. 2520 

ft 

ft 


. 03 

0. 2320 x 10 _2 ° 

1.4541 

h/5 

0. 0856 


. 06 

0. 1153 x 10~ 19 

0. 7761 

j r 

ft 

P6 

. 09 

0. 1578 x iO -13 

0. 5867 

ft 

f f 

(r = . 04) 

. 12 

0. 9120 x IO" 10 

0. 4644 

if 

1 f 


• 18 

0. 5630 x IO -6 

0. 3419 

If 

ft 


. 24 

0. 3578 x IO -4 

0.2808 

ft 

tf 
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TABLE A4. RESULTS OF E4 


Process 

h 

Error 

Time (min) 
2. 1484 

Starters 
Computed At 

Time Loss 
(min) 


. 03 

0. 5170 x 10~ 15 



. 06 

0. 2778 x 10 -13 

1. 0769 

i 


Pi 

. 09 

0. 3268 x 10" 12 

0. 7203 




. 12 

0. 1899 x 10" 11 

0. 5419 




. 18 

0. 2309 x 10" 10 

0. 3639 




. 24 

0. 1381 x 10 -9 

0. 2747 




. 03 

0. 3892 x 10“ 13 

0. 9873 

h 

0 


. 06 

0. 2265 x 10 -11 

0. 5025 

ft 

0 

P2 

. 09 

0. 1135 x 10 _1 ° 

0.4653 

n 

0 

(r = . 08) 

. 12 

0. 2672 x 10 -1 ° 

0.3067 

M 

0 


. 18 

0. 1134 x 10 -9 

0. 2100 

T J 

0 


. 24 

0 . 2 3 20 x 10 41 

0. 1394 

it 

0 


. 03 

0. 1015 x 10“ 13 

1. 1733 

h/2 

0. 0214 


, 06 

0. 1302 x 10 -11 

0. 6072 

IT 

H 

P3 

. 09 

0. 6201 x 10 -11 

0.4186 

M 

H 

' 

( r = . 04) 

. 12 

0. 7905 x 10 _1 ° 

0. 3244 

1? 

IT 


. 18 

0. 2019 x 10 -8 

0. 2303 

n 

Tf 


. 24 

0. 3918 x 10 -8 

0. 1831 

ft 

M 


TABLE A4. RESULTS OF E4 (Concluded) 


Process 

h 

Error 

Time (min) 

Starters 
Computed At 

Time Loss 
(min) 


. 03 

0. 9590 x 10~ 15 

1. 3561 

h/2 

0. 0214 


. 06 

0. 8544 x 10 -13 

0. 7914 

n 

1! 

P4 

. 09 

0. 1494 x 10" 12 

0. 5400 

M 

II 

4 

II 

o 

£ 

. 12 

0. 9297 x 10" 11 

0.3686 

IT 

II 


. 18 

0. 3043 x 10 -8 

0. 2889 

M 

II 


. 24 

0.2025 x 10" 6 

0. 2042 

IT 

IT 


. 03 ^ 

0. 8970 x 10- 15 

1. 5397 

h/5 

0. 1143 


. 06 

0. 6492 x 10~ 12 

0. 8364 

T 1 

II 


. 09 

0. 9604 x 10 -11 

0. 6020 

TT 

II 

P5 

. 12 

0. 3883 x 10 _1 ° 

0. 4397 

IT 

11 

(r = . 04) 

. 18 

0. 1866 x 10- 8 

0. 3386 

IT 

11 


.24 

' 

0.4390 x i0~ 8 

0. 3092 

IT 

IT 


. 03 

0. 1460 x 10- 15 

1. 7933 

h/5 

0. 1143 


. 06 

0. 2042 x 10" 14 

0. 9614 

IT 

II 

P6 

.09 

0. 2125 x 10 -11 

0. 6841 

IT 

11 

(r = . 04) 

. 12 

0. 7706 x 10~ 8 

0. 5458 

r r 

If 


. 18 

0. 1833 x 10~ 4 

0. 3783 

IT 

II 


.24 

0. 5610 x 10~ 3 

0. 3167 

IT 

II 


REFERENCES 


1. Collatz, L. : The Numerical Treatme nt of Differential Equations, 3rd ed. , Springer, 
Berlin, 1960. 

2. Dahlquist, Germund. : Stability and Error Bounds in the N umerical Integration of 
Ordinary Differential Equations . Almquist and Wiksells Boktrycheri AB, Uppsala, 
1958. 

3. Fehlberg, Erwin. : Numerically Stable Interpolation Formulas with Favorable Error 
Propagation for First and Second Ord er Differential Equations. NASA TN D-599, 1961. 

4. Hamming, R. W. : Numerical Methods for Scienti sts and Engineers. McGraw-Hill, 
New York, 1962. 

5. Henrici, Peter. : Discrete Variable Methods in Ordinary Different ial E quations. 

Wiley, New York, 1962. 

6. Kunz, K. S. : Numerical Analysis. McGraw-Hill, New York, 1957. 

7. Milne, W. E.: Numerical Solution of Diff erential Equations . Wiley, New York, 1953. 

8. Nordsieck, Arnold. : On Numerical Integration of O rdinary Differential Equations. 
Mathematics of Compu tation, vol. 16, no. 77, January 1962. 

9. Shanks, E. B. : Formulas for Obtaining So lutions of Differential Equations by_the 
Evaluation of Functions. Vanderbilt University, Nashville, Tennessee, 1963. 


40 


NASA -Langley, 1965 M509 



"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Adjninistration 
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